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We present a detailed account of quantum state estimation by joint maximization of the likelihood 
and the entropy. After establishing the algorithms for both perfect and imperfect measurements, 
we apply the procedure to data from simulated and actual experiments. We demonstrate that the 
realistic situation of incomplete data from imperfect measurements can be handled successfully. 
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I. INTRODUCTION 

Quantum state preparation is the first important step 
for any protocol that makes use of quantum resources. 
Examples of such protocols are quantum state teleporta- 
tion and quantum key distribution which require entan- 
gled quantum states. In order to verify the integrity of 
the quantum state of the source prepared, one carries out 
quantum state tomography on the source. Measurements 
are performed on a collection of quantum systems (elec- 
trons, photons, etc.) that are emitted from the source, 
that is, a quorum. Then, the quantum state of the source 
is inferred from the measurement data obtained from this 
ensemble. The measurements are generically described 
by a set of positive operators IX,- that compose a prob- 
ability operator measurement (POM). The procedure of 
state inference, which shall be our main focus in this arti- 
cle, is also known as quantum state estimation. If the size 
of the ensemble is infinite, the estimation procedure will 
yield the unique true quantum state of the source; this is 
the frequentist's definition of the true state, which we ac- 
cept as the best description of what the source prepares. 
However, such an ensemble is never achievable in any lab- 
oratory setting, as one can only perform measurements 
on a finite ensemble of quantum systems. As a result, the 
state estimator obtained will be different from the true 
state and depends on the details of the estimation proce- 
dure. To make statistical predictions, the corresponding 
operator p describing this estimator must be a statistical 
operator, which is positive. This will ensure that the es- 
timated probability pj = ti{pHj} for an outcome IX, of 
any set of POM is positive. We shall denote all estimated 
quantities with a "hat" symbol. 

There are two popular methods for quantum state es- 
timation: Bayesian and maximum-likelihood (ML). The 
Bayesian state estimation method [ll-Q constructs a state 
estimator from an integral average over all possible quan- 
tum states. The likelihood functional, which yields the 
likelihood of obtaining a particular sequence of measure- 
ment detection with a given quantum state, serves as 
a weight for the average. This approach includes all 
the neighboring states near the maximum of the likeli- 



hood functional as possible guesses for the unknown ptrue- 
These neighboring states are given especially significant 
weight when N is small, in which case the likelihood func- 
tional is only broadly peaked at the maximum. How- 
ever, the integral average unavoidably depends on how 
one measures volumes in the state space, and there is no 
universal and unambiguous method for that. The ML 
approach 0-05 on the other hand, simply chooses the 
estimator as the statistical operator that maximizes the 
likelihood functional. Rather than identifying a unique 
estimator, as the Bayesian approach always does, the ML 
method may only yield a convex set of estimators if the 
estimated probabilities pj are consistent with more than 
one statistical operator. If the ML estimator is unique, 
and the quorum sufficiently large, both approaches give 
the same estimator since the likelihood functional peaks 
very strongly at the maximum. 

When the measurement outcomes form an informa- 
tional^ complete set, the measurement data obtained 
will contain maximal information about the source. 
Thus, a unique state estimator can be inferred with ML. 
Unfortunately, in tomography experiments performed on 
complex quantum systems with many degrees of free- 
dom, it is not possible to implement such an informa- 
tionally complete set of measurement outcomes. As a 
result, some information about the source will be miss- 
ing and its quantum state cannot be completely char- 
acterized. The ML estimator obtained from these in- 
formationally incomplete data is no longer unique and 
there will in general be infinitely many other ML esti- 
mators which are consistent with the data. In Ref. [H, 
we briefly reported an iterative algorithm (MLME) to 
estimate unknown quantum states from incomplete mea- 
surement data by maximizing the likelihood and von 
Neumann entropy functionals. In that Letter, we as- 
sumed that the measurement detections are perfect with 
no detection losses, i.e. . = 1. The application of 
this algorithm was illustrated with examples of homo- 
dyne tomography and we concluded that, together with 
a more objective Hilbert space truncation, this approach 
can serve as a reliable and statistically meaningful quan- 
tum state estimation with incomplete data. 

In this article, we will present more details on the re- 
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cently proposed MLME algorithm and apply it to vari- 
ous other situations. First, we give a brief review of the 
mathematical formalism for quantum state estimation in 
Sec. [Ill to set the stage for the subsequent discussions. 
Next, we derive the numerical MLME algorithms respec- 
tively for both perfect and imperfect measurement detec- 
tions in Sec. IIII1 with the latter being particularly useful 
for actual experiments. We illustrate applications of the 
two algorithms with two examples in Sec. IIVI and finally 
conclude in Sec. [Vj 



II. FORMALISM OF QUANTUM STATE 
ESTIMATION 

In a tomography experiment, an ensemble of N copies 
of quantum systems, identically prepared, is measured 
using a POM which consists of positive measurement out- 
comes Ilj. For simplicity, we first assume that all mea- 
surement detections are perfect and hence . Hj = 1. 
The problem of imperfect detections will be dealt with 
in Sec. IIII Bl For each outcome, its number of occurrences 
is denoted by rij such that J2j n j = N. The likelihood 
functional £({nj};p), for a particular sequence of inde- 
pendent detections, is then 

£({n j };p) = l[p?. (1) 
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As a consequence of perfect measurement detections, 
^2 j Pj = 1. The ML procedure searches for the estimator 
Pml which maximizes C({nj}; p). For a .D-dimensional 
Hilbert space, when a POM comprises D 2 or more mea- 
surement outcomes, of which D 2 of them are linearly 
independent, it is informationally complete. In this case, 
there exists a unique estimator pML for a given set of 
measurement data {nj}. One can also define the out- 
come frequencies fj = rij/N out of these measurement 
data such that J2j fj = l- The corresponding functional 
£({rij}; p) due to this informationally complete POM will 
peak at the unique global maximum pml over the space 
of p, whereby pML is solely determined by the frequen- 
cies fj and does not depend on the total number N of 
measured copies. 

The situation is different when the POM is informa- 
tionally incomplete. In this case, there will be infinitely 
many ML estimators satisfying a smaller set of linearly 
independent constraints imposed by the incomplete mea- 
surement data. These ML estimators form a convex 
set of operators which maximize the convex functional 
£({rij}; p). Geometrically, C({rij}; p) possesses a convex 
plateau structure hovering over the space of p. The task, 
now, is to select one of these estimators for future sta- 
tistical predictions. To do this, we adopt the well-known 
maximum-entropy (ME) principle advocated by Jaynes 
0| . That is, we look for the estimator with the largest 
von Neumann entropy 

S(p) = -tr^logp} (2) 



among the convex set of ML estimators. This supplemen- 
tary step introduces a small and smooth convex hill over 
the plateau structure so that a unique maximum can be 
obtained. The corresponding MLME estimator Pmlme 
is the least-bias estimator for the given set of incomplete 
measurement data; it can be regarded as the most con- 
servative guess of the unknown quantum state out of the 
convex set of ML estimators. 

At this point, we would like to comment on the dis- 
tinction between this MLME technique and the conven- 
tional ME technique [To|, El- The ME technique takes 
the outcome frequencies fj as bona fide estimates for the 
probabilities pj and tries to search for the positive oper- 
ator 

^ ME = f V a n H ^ 
trje*^ XjUj } 

that maximizes S(p), subjected to the probability con- 
straints which are mediated by the Lagrange multipliers 
Xj. The fundamental problem with this scheme is that 
the fjs cannot always be treated as probabilities since 
there may not be any statistical operator p for which 
fj = tr {pUj}. This is due to the statistical noise which 
is inherent in the outcome frequencies arising from mea- 
suring a finite ensemble of quantum systems. Therefore, 
in such cases, the ME technique fails as there simply is 
no positive operator which is consistent with the mea- 
surement data to begin with. The MLME algorithm, on 
the other hand, looks for the unique MLME estimator by 
confining the search within the plateau region inside the 
space of statistical operators. Thus, positivity is ensured. 
In cases where the fjS are probabilities, both the ME and 
MLME schemes yield the same estimator by construction 
since the estimated probabilities pj = fj correspond to a 
statistical operator. 



III. THE NUMERICAL ALGORITHMS 

A. Perfect measurements 

Assuming that the measurement detections are per- 
fect, the likelihood functional jC({rij}; p) in Eq. (pQ) 
gives a complete statistical description of all possible 
sequences of detections for the N measured copies of 
quantum systems. Equivalent ly, one can consider the 
optimization of the normalized log-likelihood functional 
log(£({nj}; p))/N to simplify the subsequent calcula- 
tions, in view of the monotonic nature of the logarithmic 
function. The motivation for introducing the normaliza- 
tion will become clear soon. The MLME scheme can 
then be perceived as a standard constrained optimiza- 
tion problem: maximize log(£({nj}; p))/N subjected to 
the constraint that S{p) takes the maximal value S'max- 
This is equivalent to maximizing S(p) with the constraint 
that \og(C({nj}; p))/N is maximal, as discussed above. 
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The Lagrange functional for this optimization problem is 
defined as 



l(X;p) = \(S(p)-S n 



1 

N 



log £({71,}; p), (4) 



where A is the Lagrange multiplier corresponding to the 
constraint for S(p). We denote the estimator that max- 
imizes Z(A; p) by pi^\. Incidently, the functional Z(A; p) 
is a sum of two different types of entropy, up to an irrel- 
evant additive constant Y^j fj l°g fj : the von Neumann 
entropy S(p) that quantifies the "lack of information", 
and the negative of the relative entropy S({fj}\{pj}) = 
J2j fj log(fjfpj) that quantifies the "gain of information" 
from the measurement data. The scheme can now be in- 
terpreted as a simultaneous optimization of two comple- 
mentary aspects of information, with an appropriately 
assigned constant relative weight A. In addition, the nor- 
malization of \og C({nj}; p) renders the optimal value of 
A to be independent of N. 

When A = 0, we recover the Lagrange functional for 
the log-likelihood functional alone. Owing to the infor- 
mational incompleteness of the measurement data, there 
exists a convex plateau structure for the log-likelihood 
functional. As A — >> oo, the von Neumann entropy be- 
comes increasingly more significant and the resulting es- 
timator pi^A^oo approaches the maximally- mixed state 
1/D. Naturally, when A takes on a very small posi- 
tive value, the contribution from XS(p) becomes much 
smaller than log(jC({rij}; p))/N and the effect of the von 
Neumann entropy functional is only significant over the 
plateau region in which the likelihood is maximal. Fig- 
ure [D illustrates all the aforementioned points. This 
means that, in general, A should be chosen so small that 
S (pi,x) is very close to the minimum, and below which 
there are only very slight changes in the two entropy 
functionals [8[. 

Let us derive the iterative algorithm for maximizing 
X(X —> 0; p) with respect to p. After varying X(X — >• 0; p), 
we have 

8X(X^0;p) = -Atr{6plogp} + J]— 8 Pj . (5) 

j J 

The variations &pj, or 6p, have to be such that p stays 
positive after these variations. To choose their appro- 
priate forms, we first parameterize the positive operator 
p = A^/trjA^} with an auxiliary complex operator 
A. Under this parametrization, 



bp. 



bA^A + A*6A- p tr{6A*A + A^bA} 



tr{A^A} 

Substituting bp in Eq. (j6| into Eq. (j5), we have 



6J(A -)■ 0;p) = trj 
where 



tr{.AU}- 



tr{A^A} J ' 



V\ = R-l- A(logp- tr{plogp}) 



(6) 



(8) 





(a) A = 



(b) A 




(c) A = 10- 3 

FIG. 1. Schematic diagrams of X(A, p) on the space of sta- 
tistical operators. The maximally-mixed state resides at the 
center of the square base which represents the Hilbert space. 
At the extremal points of A, X(A = 0; p) — log(£({rij}; p))/N, 
with a convex plateau at the maximal value, and X(A — > 
oo; p) = XS(p). Plot (c) shows the functional with an appro- 
priate choice of value for A for MLME. An additional hill-like 
structure resulting from S(p) is introduced over the plateau, 
so that the estimator with the largest entropy can be selected 
from the convex set of ML estimators within the plateau. 



with 
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E- n i- 



(9) 



When X(X —> 0; p) is maximal, we have bX(X —> 0; p) 
and the extremal equations 



p9r = 9fy = 



(10) 



are satisfied. Therefore, to solve these extremal equations 
numerically, we iterate the equation 



Pk+i 



(4 + 64) (A k + bA k ) 
tr{(4 + 54) (A k + bA k )} 



(11) 



starting from some statistical operator pi, until k = k' 
such that the norm of pk'^k' is less than some pre-chosen 
value. We then take Pmlme = Pk' as the MLME estima- 
tor. Maximizing X(X —> 0; p) will require bX(X — ^ 0; p) to 
be positive whenever X(X — » 0; p) is less than the maxi- 
mal value. A straightforward way to enforce posit ivity is 
to set 

bA k = (§4) f = oc e^^JT ' (12) 



with e being a small positive constant. This is the 
steepest- ascent method. We have thus established a nu- 
merical MLME scheme as a set of iterative equations (|TTj) 
and (fT2|) to search for the MLME estimator using the 
measurement data obtained from perfect measurement 
detections. More compactly, the relevant iterative equa- 
tions are 



(l + e£R fc )P*(l + d*fc) 



tr{(l + dK fc )p fc (l 
Rk - 1 - A (log p k 



tr{p k \ogp k }) . 



(13) 



We note that a more efficient algorithm, using the 
conjugate-gradient method, can be derived from this 
steepest-ascent algorithm, which is the subject of a sep- 
arate discussion. 



B. Imperfect measurements 

In actual experiments, the measurement detections will 
usually be imperfect in the sense that the detection effi- 
ciency r]j of a particular measurement outcome Uj is less 
than unity. In this case, the overall outcome probabilities 



Pj = VjPj 



(14) 



will not sum to unity. Hence, we have a set of POM 
with outcomes II j = rjjUj such that G = < 1- 

A consequence of this is that the true total number M 
of copies received is not known, since only N < M are 
detected (N = M when all rjj = 1 as in Sec. IIIIAj) . 

The likelihood functional that accounts for all M 
copies of quantum systems in an experiment with im- 
perfect detections is given by 



C({nj};p) 



Ml 



N\ (M-N)l 



n 



Pi 



(l-Tj) 



M-N 



(15) 

where n = . pj < 1 . The additional combinatorial pref- 
actor arises from the indistinguishability in the ordering 
of the detection sequence resulted from losses. With the 
help of Stirling's approximation for the factorials, the 
variation of the corresponding log-likelihood functional 
is given by 



61og£({n j };p) = tr<^ NR 



M-N 



1 



5Mlog 



(1 



V 
-N 



G bp 



where R = Y^jfj^j/Pj- Adopting the concept of 
maximum-likelihood, we derive an expression for M such 
that log C({rij}; p) is maximized for any given p. This im- 
plies that the coefficient of the arbitrary bM must vanish 
and we have M = N/rj as the most-likely value of M . 



With this, the expression for jC({nj}; p) reduces to the 
simple form 



AW;p) = II? 
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(16) 



up to an irrelevant multiplicative factor, with its corre- 
sponding logarithmic variation 



61og£({n i };p) = Ntvl R 



(17) 



The additional term —bpG/n in the argument of the trace 
accounts for copies that have escaped detection. 

Defining X(X — >• 0; p) for the new POM and its 
£({rij}', p) in Eq. ([T6]) . one can derive the iterative equa- 
tions 



Pk+l 



(l + dH fc ) p k (l + e^Hfc) 

tr{(l + e^H fc ) p k (l + e^ fe )} ' 
G 



V 



(k) 



X(\ogp k - tr{p k \ogp k }) , (18) 



with 

To highlight the importance of a proper treatment of 
imperfect measurement detections, we perform a simula- 
tion on 10 3 randomly generated qubit states. Figure [2] 
compares the performance of the MLME algorithm de- 
rived in Sec. IIII Al with which we search for the MLME 
estimator by assuming that the measured data {nj} are 
all we have while ignoring the possible missing data, with 
that of the MLME algorithm derived in this section. The 
trace-class distance 



Ar = -tr{|pMLME — Ptruel} 



(19) 



is used as the figure of merit to quantify the distance be- 
tween Pmlme and ptme- The lesson here is that if one 
neglects the consequence of imperfect measurements in 
performing state reconstruction, the quality of the result- 
ing reconstructed state estimator will typically be much 
lower than that obtained from a scheme which accounts 
for this imperfection. 



IV. APPLICATIONS 
A. Time-multiplexed detection tomography 

First, we apply the MLME technique to simulation ex- 
periments on time-multiplexed detection (TMD) tomog- 
raphy [12|. For experiments of this type, photon pulses, 
of a particular quantum state, containing more than one 
photon are sent through a series of beam splitters [l3| , 
each associated with a certain transmission probability. 
Behind each of the output ports of such a series is a 
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FIG. 2. A comparison of two different schemes with 10 3 ran- 
dom qubit true states distributed uniformly with respect to 
the Hilbert-Schmidt measure. Fifty experiments were simu- 
lated for every true state, with N = 5000 for each experi- 
ment, and the respective average trace-class distances T>1^ 
were computed. The entire simulation was done with a set of 
randomly generated, informationally incomplete POM con- 
sisting of two imperfect measurement outcomes. The plot 
markers denoted by "+" represent reconstructed states using 
the algorithm in Eq. (|13|) while ignoring the imperfection of 
the measurements, and those denoted by "□" represent the 
reconstructed states using the algorithm in Eq. JT5J) that ac- 
counts for this imperfection. The significant improvement in 
tomographic efficiency with the latter algorithm is a strong in- 
dication of the importance of a proper treatment of imperfect 
measurements. 



a measurement of these outcomes only gives information 
about the diagonal entries of the statistical operator of 
the true state in the Fock basis. In order to obtain in- 
formation about the off-diagonal entries, one can, for in- 
stance, displace the current set of 2 A/ports POM outcomes 
in phase space with some complex value a k away from 
the origin using the displacement operator 



T>{a k ) 



^a k A^-atA 



(21) 



where A is the standard photon annihilation operator. 
Then, the new set of outcomes 

Uj(a k ) = ^V(a k )Ii V\a k ) , (22) 



with M being the total number of such displaced set of 
2 A/ports outcomes, do not commute with the undisplaced 
set. These displaced outcomes are suitable for a mea- 
surement that is designed to obtain information about 
the unknown true state by sampling over multiple a k s. 
Experimentally, these displaced POM outcomes can be 
realized with unbalanced homodyne detection [l5[ . 

In the simulations, four output ports, corresponding to 
a total of 2 4 = 16 POM outcomes, are considered. Two 
different true states are selected to illustrate the results 
of MLME. The first true state is chosen to be a stationary 
state of a laser given by 
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FIG. 3. A schematic diagram representing the time- 
multiplexed setup with K + 1 output ports. The TjS are the 
respective transmission probabilities for the jth beam split- 
ter. The overall efficiency for, say, the kth port is given by 
fjk = rjkO- — T k + Y\)= 



-i t . 



single-photon detector that either registers a click from 
an incoming split photon pulse, with some detection ef- 
ficiency, or does nothing. Thus, each output port has a 
certain overall efficiency fjj which is related to the rel- 
evant transmission probabilities and detection efficiency 
(See Fig. 

As a consequence of this, the POM outcomes 



U 3 =^2\ n ) C 3n H 



(20) 



will be a mixture of Fock states, with the coefficients cj n 
related to rjj [14|. If there are N ports output ports, where 
all rjjS are different, there will be 2 A/ p° rts distinct POM 
outcomes due to the binary nature of the single-photon 



detectors. In addition, Y^=i = 1 snlce the 2 

binary sequences of detection configurations constitute 
all possible events. These POM outcomes commute and 



p ss = e"" £ |n) t- ( 



(23) 



where \i is the mean number of photons [l6| . 
For the second true state, the statistical operator 
p a t = |m(c/)) (m(c/)|, where 



Ma')) 



W) + \-a') 



2(1 



-2\a'\ 2 



(24) 



is the superposition of the coherent states \a f ) and |— a 7 ), 
is chosen. The notation |m(c/)) is used to denote the ket 
for the male "Schrodinger's cat" state. See, for example, 
Ref. [17] for a survey of the family of cat states. Statisti- 
cal operators are first reconstructed from the simulated 
data. For this reconstruction, one has to decide on the 
dimension £> su b of the truncated Hilbert space for the 
reconstructions. This procedure, also commonly known 
as state-space truncation, depends on the prior informa- 
tion about the unknown state. In our case, suppose one 
knows that the mean number of photons of the source 
is fi « 4, which is the value assigned in the simulation. 
Then, one may anticipate that all the relevant informa- 
tion about the true state should be contained in a Hilbert 
space of a dimension which is close to fi. In fact, it is a 
common practice to choose £> S ub, compatible with this 
information, such that the displaced operators form an 
informationally complete POM. Then, the standard ML 
method can be applied to state estimation. We shall 
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compare the result of this approach with another, per- 
haps more objective, methodology in which we select a 
larger subspace compatible with this prior information 
and estimate the state with MLME. 

After obtaining the reconstructed statistical operators, 
the Wigner functions W(x,p) of the dimensionless posi- 
tion and momentum quadrature values, x and p respec- 
tively, are calculated in accordance with 



CO CO 



W(x,p) = 2e-\ a \ X)EHpI»> 

m—0 n—0 

2J>J <\ X + ■ 1 s S n(n- m ) p) \ m -n\ L (\m-n\) ^ ^ 



i-iy 



2J<j>! 



(25) 



where a = x + ip and 1$ (y) is the degree- n associated 
Laguerre polynomial in y of order v, for all the statistical 
operators. Here, we define j < = min{m, n} and j > = 
max{m, n}. 

To quantify the nonclassicality of the statistical opera- 
tors, we make use of the concept of nonclassicality depth 
introduced in Ref. [18|. Let us define the function 



7Z(a, r) = — / (dw) exp 

7TT 



\a/y/2- 



w\ 



P(w), 



(26) 

where w is a complex variable, (dw) denotes the integral 
measure over the real and imaginary parts of w, P(w) 
is the Glauber- Sudarshan P function, and the parameter 
r is in the range < r < 1. From the above defini- 
tion, it follows that TZ(a, r) is a continuous interpolating 
function of r from the typically singular, as well as non- 
positive, P(a/y/2) (t —> 0), to the Wigner function W(a) 
(r = 1/2), and finally to the positive Husimi Q function 
Q(a) = (a\p\a) /2tt (t -+ 1). 

The nonclassicality depth is then defined as the small- 
est value r = f, above which 7Z(a, r) > 0. Any mixture 
of coherent states is therefore a classical state since, in 
this case, f = 0. A quantum state with f > is a non- 
classical state. This measure of nonclassicality captures 
the nonclassical nature of quantum states through a one- 
parameter family of functions, which can otherwise be 
invisible to measures involving a fixed value of r, such as 
the conventional negativity of the Wigner function. Al- 
though quantifying nonclassicality with f is a somewhat 
arbitrary procedure, we adopt it here as a measure of 
nonclassicality that is not worse than other proposals. 
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(a) True state 





(b) 5-dimensional ML estimator 



(c) 11-dimensional MLME 
estimator 



FIG. 4. Density plots of the Wigner functions, in phase 
space, of various statistical operators for (a) the true state (20- 
dimensional stationary state of a laser, ji = 4) with f w 0.394, 
(b) the 5-dimensional ML estimator with f ~ 0.921 and (c) 
the 11-dimensional MLME estimator with f ~ 0.489. Here, 
brighter regions indicate the locations of larger Wigner func- 
tion values, and vice versa. The statistical operator for (b) is 
obtained using ML by assuming a 5-dimensional subspace in 
which the displaced POM outcomes are informationally com- 
plete. The statistical operator for (c) is obtained by assuming 
a larger subspace of dimension 11 using MLME. Numerous ar- 
tificial nonclassical features of the ML estimator, a signature 
of its highly oscillatory Wigner function, are manifested as 
an abnormally large value of f, an inevitable byproduct of 
state-space truncation. One can see that with MLME, extra- 
neous artifacts of the Wigner function resulted from such a 
truncation can be largely removed. 



The generalization of ([25]) to arbitrary r values, 



r 



_ |q=I oo co 

K(x,p,r) = ( m \p\ n ) 

m =0 n=0 



-(|m-n|) 




. + i sgn(n-m) \ 



\/2(l - r) 



2t(1 - t) 



(27) 



is useful for the numerical computation of f. For the 
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(a) True state (b) 8-dimensional ML estimator 




(c) 10-dimensional MLME (d) 15-dimensional MLME 

estimator estimator 



FIG. 5. Density plots of the Wigner functions, in phase space, 
of various statistical operators for (a) the true state (p a /, 
a' — 5), (b) the 8-dimensional ML estimator, (c) the 10- 
dimensional and (d) 15-dimensional MLME estimators. In 
this case, the Wigner function of the ML estimator differs 
greatly from that of the true state, an example of misleading 
information obtained via state-space truncation. A transition 
in the structure of the Wigner function occurs at D su h — 10, 
with the MLME estimator for D su h = 15 giving a more ac- 
curate estimated picture of the Wigner function of the true 
state. 

stationary state in Eq. ([23]) . Eq. ([27]) simplifies to 



K ss (x,p,r) 




(28) 

The performances of both MLME and the standard 
ML method on the true states defined in Eqs. ([23]) and 
([24]) are illustrated by the Wigner function plots of the 
respective statistical operators obtained from both meth- 
ods. These are shown in Figs. |4] and [5] The respective 
nonclassicality depths are also computed for Fig. HI For 
the state p a > , all the corresponding reconstructed statis- 
tical operators are highly nonclassical, with f = 1 [l9| 
for all of them. Hence, rather than compare the f val- 
ues, the structure of the Wigner functions for various 
reconstruction subspaces will be briefly analyzed instead 
in Fig. 
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FIG. 6. Schematic diagram of the diffraction patterns of an 
incoming light beam that is obtained from a SH wave front 
sensor. The light beam is transformed by an array of mi- 
crolenses (apertures). A CCD camera is placed at the rear 
focal plane of the array. The measurement data consist of 
the measured intensities of the beam. The intensity at the 
jth pixel, located at position Xj, behind the kth microlens 
aperture is denoted by Ik(xj). 



B. Light-beam tomography 

Finally, we make use of the MLME algorithm to recon- 
struct states of classical light beams that are measured 
using the Shack-Hartmann (SH) wave front sensor. An 
incoming light beam is transformed by a regular array of 
microlens apertures and detected in its rear focal plane 
by a charge-coupled device (CCD) camera (see Fig. [6]). 
A plane wave traversing in the transverse plane of the 
SH sensor gives rise to a detection, where the individ- 
ual diffraction patterns are centered at the corresponding 
optical centers of the microlenses. For a distorted wave 
front, the observed diffraction pattern behind the fcth mi- 
crolens aperture will be deflected by an angle Ok- Since 
the set of angles Ok is related to the local wave front tilts 
with respect to the transverse plane of the SH sensor, 
the shape of the wave front can be inferred. Clearly, this 
standard technique of wave front reconstruction fails in 
the presence of imperfect coherence, where the notions 
of "wave front" and "optical phase" are no longer well- 
defined and a more general description of the state of the 
light beam is necessary. 

Recently, an alternative theory for SH detection, based 
on the principles of quantum state tomography, has been 
introduced. It was shown that a complete characteriza- 
tion of a beam of light is possible from the measurement 
data obtained with the SH sensor under certain assump- 
tions with regard to the aperture profiles [20]. Analo- 
gously to quantum states, we can describe a coherent 
beam (mode), with a complex amplitude ip(x), by a ket 
|^), such that ip(x) = (x\ip). It should be understood, 
that this ip(x) is not a quantum mechanical probability 
amplitude, but a mathematical symbol with analogous 



properties that we exploit. At the focal plane of the fcth 
microlens aperture, the amplitude ^' k {x) of the trans- 
formed beam is given by 



J dx' hk(x 



x')a k (x')^(x'), (29) 



where a k {x) is the aperture function of the fcth microlens 
aperture and the response function hk(x) describes the 
free propagation from the fcth microlens to the SH sensor. 

Now, suppose a generic partially coherent beam is de- 
tected by the SH sensor. We can describe the state of 
such a beam with a coherence operator p co h- When us- 
ing a computational basis of orthonormal modes \ip n ), we 
have 



Pcoh = ^2\^m) Pmni^n 
mn 

By defining the aperture operator 



(a) 



dx' \x') CLk{x') (x f 



(30) 



(31) 



for the fcth microlens aperture and the unitary propa- 
gation operator where (x\ U k \x') = h k (x — x f ), that 
describes the free propagation from the fcth microlens to 
the SH sensor, the representation of the corresponding 
transformed state p' coh . 



Pcoh 



U k Ml a) p coh Ml a) ut 



>P™(ll>n\Ml a) Ut 



(32) 



on the focal plane of the apertures follows from the lin- 
earity of optics transformations. The intensity Ik(xj) at 
position Xj [2l| on the rear focal plane of the fcth aperture 
is 



I k {xj) = {xj\p f coh \xj) 



^ mn 



(33) 



where ip f n k (xj) = (xj\ip f n k ) are the complex amplitudes 
of the transformed light beam obtained from the ampli- 
tudes tjj n (xj) = (xj\ip n ) of Eq. ([29]) . Since p co h possesses 
all the properties of a statistical operator, the MLME 
technique can be used to estimate the true coherence op- 
erator Pcoh 6 °f a partially coherent beam. To this end, 
we need to compute the corresponding POM describing 
the measurement outcomes of the SH sensor. By relating 
I k (xj ) to the corresponding probabilities of the outcomes 

U k(Xj)=J2mn Wm)^nm{Xj) «|, We have 



h{xj) = tr{p coh IL k (xj)} 

— Pmn ^-k,nm(Xj) . 



(34) 




FIG. 7. Experimental set-up involving a 
(SMF), a spatial light modulator (SLM), 
(A) and a Shack- Hart mann (SH) sensor. 



single-mode fiber 
an aperture stop 



Comparing Eqs. ([33]) and (fM]h the positive operator de- 
scribing the detection outcome at the jth pixel of the 
CCD camera behind the fcth aperture is given by 



(35) 



As an illustrative example, the POM outcomes con- 
sidered in this section are commuting operators in the 
infinite-dimensional Hilbert space with regard to the co- 
herence operators. Equivalent ly, the aperture functions 
for the respective microlenses do not overlap in position. 
This is a special case of a more general theory on Shack 
Hartmann detection, which will be discussed at length in 
another upcoming article. 

In the experiment, a controlled preparation of optical 
beams is realized using the principles of digital hologra- 
phy [22|. Figure [7] shows the set-up. The essence of the 
beam preparation lies in the numerical construction of a 
digital hologram that is programmed to produce a super- 
position of a reference plane wave and a beam with the 
true state of interest. This is achieved with the help 
of an amplitude spatial light modulator (OPTO SLM) 
with a resolution of 1024x768 pixels. The hologram is 
then illuminated by the reference plane wave that is con- 
sidered in the superposition. To approximately produce 
this plane wave, a collimated Gaussian beam is generated 
by placing the output of a single-mode fiber at the focal 
plane of a collimating lens. In this way, the digital holo- 
gram can be fully situated at the center of the collimated 
Gaussian beam of a larger beam waist, where this beam 
can then be approximated to be a plane wave with high 
accuracy. The resulting diffraction spectrum, after illu- 
minating the digital hologram with the collimated Gaus- 
sian beam, involves several diffraction orders, of which 
only one contains useful information about p^f. To fil- 
ter out the unwanted diffraction orders, a 4-/ optical 
processor, with a small circular aperture stop placed at 
the rear focal plane of the second lens, is used for this 
purpose (the aperture stop in Fig. [7]). The resulting light 
beam with the state pj;™ is then focussed at the rear focal 
plane of the third lens. This completes the preparation 
stage. 

The measurement of the light beam involves a Flex- 
ible Optical SH sensor with 128 microlenses that form 
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FIG. 8. CCD image for the state pl^. The relevant part of 
the SH readout used for the beam reconstruction is shown. 
Contributions from the individual SH apertures are indicated 
by bright spots, with each spot made up of multiple pixels. 
Note that the two void regions correspond to the phase sin- 
gularities of the state p™*. This hints that pJX « p™%. 




a hexagonal array. Each microlens has a focal length 
of 17.9mm and a hexagonal aperture with a diameter of 
0.3mm. The signal at the focal plane of the array is de- 
tected by a uEye CCD camera that has a resolution of 
640x480 pixels, with each pixel being 9.9|Ltmx9.9|LLm in 
size. 

The aforementioned set-up is used for generating 
and analyzing low-order Laguerre-Gaussian (LG) modes. 
The LG modes can serve as important resources in quan- 
tum information processing [23[. In this experiment, only 
LG modes with no radial nodes are considered. Such 
modes form a one-parameter orthonormal basis, where 
the modes are specified by the orbital angular momen- 
tum quantum number I. In polar coordinates, the rele- 
vant part of the complex amplitude of a LG mode, for a 
fixed Z, is given by 



FIG. 9. MLME state estimation from informationally in- 
complete data for D su h = 9. The real (left) and imaginary 
(right) parts of the reconstructed coherence operator p^ ME 
are shown. The reconstruction subspace is spanned by the 
modes LGj, with / = 0, 1, . . . , 8. In this case, 56 out of 91 in- 
dependent outcomes, required for complete characterization 
of ^oh e 5 are not accessible, yet the MLME estimator p^ ME 
is close to p co ^ 5 with a fidelity of 92%. 



cases where state reconstruction on informationally com- 
plete subspaces gives unsatisfactory results, the MLME 
approach can be used on the informationally incomplete 
data to give reasonable estimators on a larger subspace, 
as illustrated in Fig. [9] 



(s,<p\LGi) ocs'e^e 



(36) 



Nonzero values of I give rise to helical wave fronts, for 
which each photon carries an orbital angular momentum 

of in. 

For the source of light beams, we would like to prepare 



the state 



sup 
^coh 



sup 1 1 



where 



|Vsup) = (|LG )-|LGi)i-|LG 2 )) 



V3' 



(37) 



using the OPTO SLM. In the presence of experimental 
imperfections, however, the true state pj;™ prepared this 
way will not be exactly the same as p^ . After measuring 
this beam with the SH sensor, the data are processed 
using the MLME algorithm in Eq. (fT8|) to obtain the 
estimator p^ ME for p%™, since G < 1. To quantify 
the quality of p^ ME , we investigate the fidelity between 
and pZl 

Figure [8] shows the CCD image for the state p^^. Each 
aperture gives rise to a bright spot in the CCD image. To 
maximize the signal-to-noise ratio, only the pixel with the 
highest intensity within each spot is selected as a mea- 
surement datum. The set of intensities, corresponding to 
maximum-intensity pixels, constitute the measurement 
data to be used for state reconstruction. In our case, 
the corresponding POM consists of 35 linearly indepen- 
dent outcomes described by Eq. ([35]) . This measurement 
is, therefore, informationally complete for Z} su b < 5. In 




FIG. 10. Average fidelities, computed over 50 random choices 
of computational bases, of the estimators for different dimen- 
sions D su h of the reconstruction subspace. The unfilled (filled) 
circular plot markers correspond to informationally complete 
(incomplete) tomography, respectively. 

So far, the procedure of state-space truncation is per- 
formed in the basis of the LGi modes. In this basis, when 
Peon 6 i s known to be quite close to /^oh> ^ ne truncation of 
modes of higher orders will not result in a great loss of 
reconstruction information, as implied by the structure 
of p s ™^ in Eq. ([37]) . The situation will be very different 
when there is no such prior knowledge about pl^, except 
for the fact that the possible values of I lie in a certain 
range. In this situation, there is no appropriate strategy 
to choose a computational basis in which the state-space 
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truncation can be done effectively and justifiably. More 
generally, estimating the unknown state p 1 ^ on a trun- 
cated subspace will, as a rule, result in missing important 
reconstruction information and this will lead to strongly 
biased estimators. A remedy for this problem is to per- 
form state reconstruction on a sufficiently large subspace 
that is compatible with the knowledge about the range 
of values of I. 

To emphasize this point, we simulate the following sce- 
nario: 

• The set of measurement data, obtained from the 
CCD image shown in Fig. [8j is distributed to 50 
parties. The possible values of I for the true state 

are known to lie in the range I G [0, 7]. 

• Each party selects a computational basis and esti- 
mates the state of the beam for Z} su b = 3, 4, . . . , 8 
using either the ML (for D suh < 5) or the MLME 
algorithm (for Amb > 5). 

• The reconstructed estimators for the six values of 
D su b are reported by each party and the average 
fidelity of the estimators for every value of -D su b 
are calculated. 

A typical outcome of this scenario is shown in Fig. [TUJ As 
can be seen, performing state-space truncations in order 
to reconstruct p 1 ^ with an informationally complete set 
of data generally leads to low fidelities in the estimators. 
Increasing the number of degrees of freedom and using 
the MLME algorithm to cope with the completeness issue 
seems to be a much better strategy. 



V. CONCLUSION 

We derived the iterative algorithms for information- 
ally incomplete quantum state estimation respectively for 
perfect and imperfect measurements. Next, we applied 
these algorithms to time-multiplexed detection tomogra- 
phy and light-beam tomography. From these two appli- 
cations, we learned that one should better not restrict 
the state reconstruction to a subspace in which the rele- 
vant measurements are informationally complete. Doing 
so can result in reconstruction artifacts that originate 
in the state-space truncation and may result in inaccu- 
rate estimators for the unknown true state. Instead, one 
should perform the reconstruction on a larger subspace, 
with additional unsampled degrees of freedom, that is 
compatible with any prior information about a given un- 
known state. Such a more objective way of state estima- 
tion results in a much better tomographic quality of the 
reconstructed estimator. 
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